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We combine the coherent modified Redfield theory (CMRT) with the equation of motion-phase 
matching approach (PMA) to calculate two-dimensional photon echo spectra for photoactive molec¬ 
ular complexes with an intermediate strength of the coupling to their environment. Both techniques 
are highly efficient, yet they involve approximations at different levels. By explicitly comparing with 
the numerically exact quasi-adiabatic path integral approach, we show for the Fenna-Matthews- 
Olson complex that the CMRT describes the decay rates in the population dynamics well, but final 
stationary populations and the oscillation frequencies differ slightly. In addition, we use the com¬ 
bined CMRT-I-PMA to calculate two-dimensional photon-echo spectra for a simple dimer model. We 
find excellent agreement with the exact path integral calculations at short waiting times where the 
dynamics is still coherent. For long waiting times, differences occur due to different final stationary 
states, specifically for strong system-bath coupling. For weak to intermediate system-bath couplings, 
which is most important for natural photosynthetic complexes, the combined CMRT-I-PMA gives 
reasonable results with acceptable computational efforts. 


I. INTRODUCTION 

Photosynthesis is the process used by plants and bac¬ 
teria to convert the energy of sunlight into chemical en¬ 
ergy in order to fuel the organism’s activities. In the ini¬ 
tial steps of the photosynthetic process, pigment-protein 
complexes complete the light-energy transfer and charge 
separation with a near unity quantum efficiency. Light¬ 
harvesting molecular complexes achieve the energy trans¬ 
fer by using an array of light-harvesting pigments that 
absorb the energy to form excitons and funnel the excita¬ 
tion to the reaction center [l|. To investigate the energy 
transfer dynamics in the first steps of photosynthesis on 
the fs time scale, ultrafast spectroscopic tools are avail¬ 
able. Among the techniques used, two-dimensional (2D) 
photon echo spectroscopy is a powerful tool which allows 
for direct mapping of the excitation energy pathwws as a 
function of absorption and emission wavelength [^ . It is 
particularly useful in examining photosynthetic systems 
in which the manifold of electronic states is closely spaced 
and broadening through static disorder yields highly con¬ 
gested spectra. Recent experimental 2D electronic spec¬ 
troscopic studies of the Fenna-Matthews-Olson (FMO) 
complex observed coherent beating signals and, thus, 
raised interest in the interplay between energy trans¬ 
fer, long-lived quantum coherence in photosynthetic pro¬ 
cesses [3| , and low-frequency vibrations of the molecular 
back bone. Also in photoactive marine c^ptophyte al¬ 
gae 0 , the light-harvesting complex LH2 Q of rhodobac- 
ter sphaeroides, and in the reaction center SH of the 
Photosystem II, long-lived oscillations have been experi¬ 
mentally observed at low (77 K) and room temperatures 
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(300 K), indicating a strong vibronic coupling in these 
systems. 

To analyze the experimental findings in such large and 
complex photoactive molecular complexes, a thorough 
comparison with theoretical calculations is essential, in 
order to arrive at a reliable interpretation of the mea¬ 
sured 2D spectra. Since it is a difficult and computa¬ 
tionally demanding task to determine 2D optical spec¬ 
tra, often only the population dynamics of exciton states 
is calculated. For the FMO complex, a rather small 
light-harvesting complex, the hierarchical equations of 
motion Q were applied and quantum oscillations were 
observed on the time scale of the 2D experiments em¬ 
ploying an environmental Debye model spectral density 
with rather small reorganization energy [9|. Employing 
the numerically exact quasiadiabatic propagator path in¬ 
tegral (QUAPI) allowed to use a more realistic measured 
environmental spectral density. However, this resulted 
in a decay of the decoherence faster than experimentally 
observed [Eini- This spectral density could, more re¬ 
cently, be employed to calculate the 2D spectra of FMO 
with the hierarchy equation [l^ and a reasonable agree¬ 
ment between theory and experiment could be achieved. 
The calculations of QUAPI and the hierarchical equa¬ 
tions of motion treated the coupling of the complex to en¬ 
vironmental fluctuations numerically exactly. However, 
the computational effort is immense, which makes the 
simulation of larger light-harvesting molecular complexes 
(which contain, typically, dozens to hundreds of excitonic 
subunits) virtually impossible. The need for a highly ef¬ 
ficient numerical tool to calculate 2D optical spectra of 
large molecular complexes with a reasonable numerical 
effort and a satisfactory accuracy still exists and it is 
expected to continue to increase. 

Given their complex molecular structures, for the cal- 
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culation of 2D spectra of large light-harvestors, approx¬ 
imate schemes are usually inavoidable. Standard Red- 
field equations , which invoke a lowest-order Born and 
a Markovian approximation, are good at weak system- 
bath coupling, but fail for strong coupling. The regime 
of intermediate system-bath coupling as present for the 
exciton dynamics in photoactive complexes [T^ is typi- 
call y a lso not properly treated within Redfield equations 
[TsIr^ . Thus, the modified Redfield theory (MRT) has 
been widely used for the description of energy transfer 
processes of large molecules (isl - l^ . In this approach, 
that contribution of the system-bath coupling Hamilto¬ 
nian, which is diagonal in the eigenbasis of the system, 
is included fully, while a second-order perturbative ap¬ 
proximation is used for the off-diagonal coupling terms. 
The equation of motion of the MRT includes a popu¬ 
lation transfer within the reduced density matrix, but 
the accompanying population-transfer induced dephas¬ 
ing is neglected. The accuracy of the MRT in view of 
the dynamics of the reduced density matrix has been an¬ 
alyzed in detail [2^. Moreover, MRT has been shown to 
have a somewhat wider range of applicability when com¬ 
pared to both the original Redheld and Forster theory 
[l9|. Also, linear absorption spectra for an ensemble of 
B850 rings have been determined which shows that MRT 
includes non-Markovian effects which clearly show up in 
the high-energy part of the static absorption lineshapes 
[ 2 ^. Different energy transfer components of LHCII 
trimer and phycoerythrin 545 have been revealed using 
MRT by simultaneous quantitative fits of the absorption, 
linear dichroism, steady-states fluorescence spectra, and 
transient absorption kinetics upon excitation at different 
wavelengths (^ . 

A more refined description of the quantum dissipative 
exciton dynamics is achieved in this work upon observing 
that the population-transfer induced electronic dephas¬ 
ing can be efficiently included in the quantum master 
equation. The off-diagonal terms in the quantum master 
equation now include the decoherence of excited states 
and electronic dephasing between ground and excited 
states by exploiting the relation l/r 2 = l/2ri + 
to estimate the different contributions to the dephasing 
rate, where T 2 is the transverse relaxation time and Ti, 
T 2 are the longitudinal relaxation time and pure dephas¬ 
ing time, respectively. While working out the details with 
the results reported in this paper, this extended quantum 
master equation has also been independently put forward 
very recently in Ref. [ 2 ^ and has been named the coher¬ 
ent modified Redfield theory (CMRT). To avoid confu¬ 
sion, we use to this nomenclature also here. 

For calculating 2D photon-echo spectra, essentially two 
different approaches are available. On the one hand, the 
response to the sequence of applied laser pulses can be 
calculated by evaluating the third-order optical response 
function [^. Modified Redfield theory was successfully 
applied to simulate the 2D spectra of the double-ring LH2 
aggregate of purple bacteria including both the B800 and 
the B850 ring [30j . Comparing experimentally measured 


and theoretically calculated results of 2D spectra revealed 
that excitation energy transfer through the LHCII hap¬ 
pens on three time scales: sub-lOOfs relaxation through 
spatially overlapping states, several hundred femtosecond 
transfer between nearby chlorophylls, and picosecond en¬ 
ergy transfer steps between layers of pigments [3l[ . More 
recently, 2D spectra of the reaction center in the Pho¬ 
tosystem H were calculated with MRT at low tempera¬ 
ture, and the charge separation process was investigated 
by MRT including charge transfer states in the model 

An alternative approach to calculate 2D optical spec¬ 
tra, which is especially useful when finite durations of 
the laser pulses as well as pulse overlap effects are taken 
into account, is the equation of motion-phase matching 
approach (PMA) [s^. Using the PMA in combination 
with the conventional Redfield equations, 2D spectra of 
a single FMO subunit were studied, and the signature 
of energy transfer was revealed by well-resolved peaks in 
the simulation with adjustable pure dephasing of exciton 
states [ 13 . 

Although MRT is used to tackle many different prob¬ 
lems in the study energy transport in photosynthetic 
complexes, no investigation of its reliability in calculat¬ 
ing nonlinear and, specifically, 2D optical spectra is at 
hand. In this work, we first verify the CMRT approach 
by comparing the population dynamics of FMO exciton 
states with numerically exact results of the QUAPI ap¬ 
proach. In addition, we combine the CMRT with the 
PMA to calculate 2D photon-echo spectra for a sim¬ 
ple dimer model. Again, the results of CMRT-l-PMA 
are benchmarked against numerically exact results of the 
QUAPI approach. For the long-time steady state dy¬ 
namics, the CMRT-l-PMA and QUAPI simulations show 
differences for intermediate and strong system-bath cou¬ 
pling. However, for intermediate coupling, as it is typical 
in photosynthetic complexes, the short time dynamics in¬ 
cluding dephasing times and coherent beating frequencies 
are well described by CMRT-l-PMA. Hence, an efficient 
numerical scheme to calculate 2D photon-echo spectra 
with a reasonable computational effort is available. 

The remainder of this paper is organized in the fol¬ 
lowing way. In the next Section |TT1 we briefly intro¬ 
duce the model of the FMO complex, for we compare 
the performance of the the CMRT-l-PMA in calculat¬ 
ing a non-trivial population dynamics. Additionally, the 
dimer model is introduced for which we compare results 
of 2D photon echo spectra. A brief description of the 
CMRT-l-PMA and QUAPI for calculating the reduced 
density matrix and 2D spectra is given in Section IHIl In 
Section HVl we give the results of the comparison, and a 
thorough discussion is appended in Section |Vl before we 
finish with a conclusion. 
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II. MODEL 


In the framework of open quantum systems the Hamil¬ 
tonian of the complete system H can be decomposed as 
four parts 


H = Hs + Hsb + Hb + Hrem 

N N 

^ ^ “f ^ ^ ^ ^ “f ^71^772)5 

m—1 m—1 n<m ( 

N 2 . 

Hb = J2 E(^ + 

m—1 j—1 


where is the on-site transition energy and Jnm is 
the intermolecular coupling. JV is the total number of 
monomers. JV^ is the number of bath modes coupled to 
molecule to, which we will take to be infinity. Xmj and 
Pmj are the mass weighted position and momentum of the 
jth harmonic oscillator bath mode with frequency uJmj- 
The interaction term Hsb = induces the 

coupling between system and bath. It is assumed to be 
separable such that K^. only acts on the system sub¬ 
space and ^m(x) only on the bath degrees of freedom. 
In the following we further assume a linear relation be¬ 
tween bath coordinates and the system. The system-bath 
interaction is then given as 

J^SB = ^ E Cfjij^mj 1 ^2^ 

m j 


and we furthermore restrict our considerations to pure 
electronic dephasing only, i.e. Km = The renor¬ 

malization term is 
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This term compensates for artificial shifts of the system 
frequencies due to the system-bath interaction. 

The influence of the bath is fully described by its bath 
spectral density 


Tm(w) = 7 ry^ 
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‘^^mj^mj 
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with reorganization energy A and high-frequency cut off 
7 . We assume the bath at each monomer to be indepen¬ 
dent (no cross-correlation between baths [s^) but with 
identical spectra, i.e. Jni^j) = Jml^)- 

Laser pulses acting on the exciton complex result in the 
addition of a system-field interaction, i.e. Hs —t Hs + 
F{t), which is defined within the dipole approximation 
according to 

F{t) = -XE{t)+H.c., (5) 

with the electric field of the laser pulse E{t) and the elec¬ 
tronic transition dipole operator fi of the exciton system 

N 

= X + X'^ with X = ^ 

m—1 


pLm determines the dipole strength and direction of the 
TOth monomer. 


A. Dimer 

The dimer is modeled by two monomers with site ener¬ 
gies Cl = —50 cm“^ = —€2 and a coupling J = 150 cm“^. 
The two dipole moments are considered to be perpendic¬ 
ular to each other, i.e. -L /.t 2 - 4^or the bath spectral 
density as specified in Eq. we choose A = 50 cm ^ 
and 100 cm“^, and set 7 = 100 cm“^. 


B. FMO 


We model a monomer subunit of the FMO trimer by in¬ 
cluding seven bacteriochlorophyll a molecular sites. The 
single excitation subspace is described by a Hamiltonian 


Hp 
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435 39.7 
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in units of cm“^, where we use the site energies and dipo¬ 
lar couplings determined by Adolphs and Renger for 
the FMO complex of Chlorobium tepidum. Bath param¬ 
eters are chosen as A=35 cm~^ and 7=53 cm“^ following 

Ref.[|. 


III. METHODS 


A. Coherent modified Redfield theory (CMRT) 


The CMRT can be derived from the Nakajima-Zwanzig 
equation [ 2 ^ using a scheme for the separation of the to¬ 
tal Hamiltonian which does not treat the whole system- 
bath interaction term Hsb perturbatively [l5|[37|. In¬ 
stead, the Hamiltonian is separated according to 


Ho = Hs + Hb + ^ \p) (/i| Hsb |m) (mI i 

77'= E Im) (mI 77sb k) , 


( 8 ) 


where |/i) are eigenstates of Hs and H' is the off-diagonal 
term of the system-bath interaction part in the exciton 
basis. In this basis. Ho is diagonal and the matrix ele¬ 
ments read 


(mI 77o Im) — + Hb{p), (9) 

where is the ^th excitonic level of the system Hamil¬ 
tonian and 
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is the weighted reorganization energy. Moreover, 


HbM 






(mI Kk Im) 


( 11 ) 


describes a bath of harmonic oscillators with mass 
frequency wj and momentum pj, shifted due to the cou¬ 
pling with the exciton state |/i). 

In addition to the redefinition of the system and the 
bath Hamiltonian, one has to define a different type of 
projection operator which only projects on the diagonal 
part of the system density matrix in the eigenstate basis. 
This is achieved by 

N 

P = ^P^ with P^- = i?^,tr{|/i) (^1 •}, (12) 

/i—0 


where P^ is the projector onto the ^th excitonic state and 
= exp{—/3HB{fJ-))/Z^^ is the equilibrium density ma¬ 
trix of the bath when the system is in the excitonic state 
1^). Here, = trexp(-/3HB(/i)) with P = 1/(A:_bT) 
and T being the temperature. 

Inserting these definitions into the Nakajima-Zwanzig 
equation, determining H' up to second order and invok¬ 
ing the time-dependent population transfer rate, one ob¬ 
tains an equation of motion for the population transfer 
terms in the form 

d 

(13) 

with the population transfer rates [T^ 

Rtj.tj.i^i^{t) = 2Re / dTtr{\v) {v\ exp{-iHoT)H’ \p.) {p.\ 
Jo 

X R>^gexp{iHQT)H'}, 

= 2Re/ d'rexp[-za;^,,r - g^^^^(r) - g^j,j,j,(T) (14) 

Jo 

+ {j ) + 

X {g^uy^ir) - [g ij) ~ + 2iA y^yy\ 


Here, = C/i — The lineshape function g^v^i'v'{t) 
can be written as the two-time integral of the bath cor¬ 
relation function according to 


g^u^l.'u'[t) = ^{p,\Kk\v) {p'\Kk\v') [ dr f dT'Cir'), 

° ° (15) 

withCt = / —- 7- 

7_oo Ti" - 1 


To obtain Eq. da, we have used the cumulant expansion 
up to second order in the system-bath coupling and have 


taken the independent bath model into account. The 
absorption lineshape within the CMRT is given by 


nOO 

I{uj) = Re'^dfB / dtexp[i{uj - ujfBo)t - gfj.fj.fj.fj.it) 
Jo 


-te r^/^/^-(^)]- 


(16) 
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as detailed in Ref. [^ . 

Up to this point, Eq. (fT^ constitutes the modified Red- 
field theory, as developed and applied in Refs. [Tsl - I^ . 
Based on the population transfer term in Eq. a , we ex¬ 
tend the quantum master equation by including also the 
coherence (or, off-diagonal) terms of the reduced density 
matrix. The resulting coherent modified Redfield quan¬ 
tum master equation now reads 


^p{t) = -i[H + Fit), pit)] - 5R{p(t)} , (17) 


where Fit) is the time-dependent system-field interaction 
term. 

The relaxation and dephasing operator JR{p(t)} now 
also includes diagonal and off-diagonal terms. The diag¬ 
onal part of the relaxation operator, which was desribed 
in Ref. [s^, reads 


»(/>(*)},, = E [Rfi^yy (Op uu Ruu flflij) P fjfj\ ■ 

U^fJ, 


(18) 


The off-diagonal terms 5R{p(t)}^j/ are now included in 
order to describe decoherenece of excited states and elec¬ 
tronic dephasing between the ground and excited states. 
Here, we use an efficient way to obtain the associated 
rates by exploiting the relation l/r 2 = 1/211 + I/T 2 
to estimate the different contributions to the dephas¬ 
ing rate. T 2 is the transverse relaxation time, Ti, T 2 
are the longitudinal relaxation time and pure dephasing 
time, respectively [l^- In detail, 1/Ti = 

'^e=itv and l/T^ Is given by the first derivative of 

lineshape function gfj,fj,uvit)■ Therefore, the off-diagonal 
terms of the excited states and between the ground and 
excited states can be written as 


'^{P[t)}fju = 


E! Rvueeit) j + 

e^fj. j 

gfjfjvuit)] Pfjvit), 

.1 


»Wf)).»={=(E 


m^n 

gfifififiit)}Pfioit)- 


n^m 


This extended quantum master equation has also been 
independently put forward very recently in Ref. and 
has been named the coherent modified Redfield theory 
(CMRT). It is an efficient, but approximate way to take 
into account population transfer and dephasing on the 
same footing. 
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B. QUAPI 

The quasiadiabatic propagator path integral [10, |4l| is 
a numerically exact approach to determine the influence 
of environmental fluctuations on the system dynamics 
within a open quantum systems approach. Specifically, 
QUAPI determines the time dependent reduced statisti¬ 
cal operator p{t) of the system. It is well established in 
the literature and we only briefly summarize the central 
features in the following. The algorithm is based on a 
symmetric Trotter splitting of the short-time propagator 
/C(tfc+i,tfc) for the full Hamiltonian into two parts, one 
depending on the system Hamiltonian, and one involving 
the bath and the coupling term. The short-time propa¬ 
gator determines the time evolution over a Trotter time 
slice 5t. The discrete time evolution becomes exact in 
the limit 5t —>■ 0. For any finite dt, a finite Trotter error 
occurs which has to be eliminated by choosing 5t small 
enough to achieve convergence. On the other side, the 
environmental degrees of freedom generate correlations 
which are non-local in time. For any finite temperature, 
these correlations decay on a time scale denoted as the 
memory time scale. The QUAPI scheme defines an aug¬ 
mented reduced density tensor, which lives on this full 
memory time window. Then, an iteration scheme is es¬ 
tablished in order to extract the time evolution of this 
object. All correlations are completely included over the 
finite memory time Tmem = K6t but are neglected for 
times beyond Tmem- One increases the memory param¬ 
eter K until convergence is found. The two strategies 
to achieve convergence, i.e., minimize 5t but maximize 
Anem = K5t, are naturally counter-current, but never¬ 
theless convergent results can be obtained in a wide range 
of parameters, including the cases presented in this work. 


C. Two-dimensional electronic spectroscopy 


In the equation of motion-phasing matching approach 
(PMA), the polarization in the photon-echo direction is 
calculated by simultaneously propagating three auxiliary 
density matrices pi , p 2 and pa, thereby employing also 
the rotating wave approximation . The time evolution 
equations are given by 

- Vi{t,ti) - V^{t,t2) - V^{t,t:i),pi{t)\ 
-lR{pi(t)}, (20) 

^P2(i) = -i[H - V^{t,t2),p2{t)\ - 5i{p2(f)}, 

^P3(i) = -i[H - V^{t,t3),p^{t)] - ?fi{p3{t)}. 

where Va{t,ta) = XEa{t — to) = XEa{t — and 

Ea(t — to) = exp(—41n2(t — ta)^jT^'), Tp is the pulse 
duration. To obtain the third-order 2D signal, the po¬ 
larization in the phase matching direction is evaluated 


as 

PpE{ti,t2,h) = E'^'‘"^{X{pi{t) - P2{t) - Psit))) +C.C., 

( 21 ) 

where the bracket (...) denotes the trace. Experimen¬ 
tally, in the limit of ideal detection, the heterodyne 
photo echo signal is proportional to the polarization 
^pe(U,^ 2 )^ 3 )^)) where t is the detection time. There¬ 
fore, the ideal total 2D signal can be expressed as 

/ OO poo 

dr / 

-OO j — OO (^^) 

X iPpEir, T,t), 

where r, T and t denote coherence time, population 
(waiting) time and detection time, respectively, t = 
^2 — U, T = ts — t 2 - The coherence time corresponds 
to a period in which the system is coherently evolving 
after the first interaction with the optical field. The sec¬ 
ond interaction with the field creates population states 
and the third interaction recovers the coherence again. 
The Fourier transform in Eq. (I22|) is always performed 
over the coherence time r and the detection time t. The 
corresponding frequencies Wt, are often referred to as 
absorption and emission frequencies, respectively. In ad¬ 
dition, Gaussian laser pulses have been assumed for a 
realistic detection scheme, which have the form 


3 

E{t) = + c.c. 

a—1 


(23) 


where A, ta, ka, and uj are the amplitude, envelop cen¬ 
tral time, wave vector and frequency of the pulses and 
Tp characterizes the pulse duration. Note that all the 
pulses are assumed to have the same lineshape, carrier 
frequencies and durations in this paper. 


IV. RESULTS 

A. Population dynamics of the FMO complex 

In order to verify the reliability of the CMRT, we 
present the population dynamics of the EMO complex 
calculated by CMRT and compare the results to those 
obtained by the numerically exact QUAPI method. In 
Eig. [U the population dynamics of selected EMO sites is 
shown for T = 77 K for two different initial conditions. 
In Eig. [T^), we assume the energy transfer to start from 
site 1. We monitor then the full transfer which involves 
all seven EMO sites. Eor simplicity, we only show the 
population dynamics of the sites I, 2, and 3. Alterna¬ 
tively, the energy transfer may be assumed to start from 
site 6, see Fig. There, we depict the population dy¬ 
namics of the relevant sites 3,5, and 6. We observe that 
the oscillatory behavior of the populations is captured by 
both approaches. Both also yield the same decay rates 
and periods of oscillations. However, a phase shift of the 
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FIG. 1. Population dynamics of selected FMO sites. In a), 
the population of sites 1 (black), 2 (red) and 3 (green) with 
the initial condition p(l, 1) = 1 is shown, while in (b) the 
populations of sites 3 (green), 5 (blue) and 6 (magenta) with 
the initial condition p(6, 6) = 1 is depicted (symbols: QUAPI, 
full lines: CMRT) for the parameters as given in the text. 


oscillations occurs between the CMRT and QUAPI re¬ 
sults. Energy transfer is believed to be related to the 
population of the FMO site 3 (green symbols and lines) 
which has the lowest energy in the FMO monomer. In 
our comparison, CMRT slightly overestimates the pop¬ 
ulation transfer efficiency towards site 3. All in all, the 
CMRT results for the FMO exciton population dynam¬ 
ics are in good agreement to numerically exact QUAPI 
results. Since the system-bath coupling parameters of 
the FMO complex are typical for natural photosynthetic 
units, we conclude that CMRT is a useful tool to study 
their exciton dynamics. 


B. Two-dimensional spectra of a dimer system 

To obtain 2D spectra, we combine the CMRT next 
with the PMA. This constitutes a very efficient ap¬ 
proximate numerical tool whose reliability is assessed 
by a comparison with 2D spectra obtained by QUAPI. 
Since 2D spectra involve extended numerical calcula¬ 
tions, QUAPI results are available only for small model 
systems with present day hardware technology. For such 
a comparison, we present the calculated results for the 
dimer model. It allows us to study energy transfer and 
dephasing (homogeneous broadening) as building blocks 
of the exciton dynamics in larger molecular compounds. 
It can still be treated by QUAPI with reasonable numer¬ 


ical effort. 

Fig. [5] (left) shows 2D spectra of the dimer calculated 
by CMRT-t-PMA for A=50 cm“^ and the other parame¬ 
ters as indicated above. They are compared to QUAPI 
results (right column in Fig. [5]) for waiting times T = 0 
is, 50 fs, 100 fs and 500 fs. These 2D spectra show two di¬ 
agonal peaks (labeled A, B) which correspond to the two 
exciton states. Moreover, two cross peaks (labeled C and 
D) arise due to the excitonic coupling between them. For 
the sake of simplicity and clarity of the comparison, inho¬ 
mogeneous broadening and the rotational averaging for 
different laser polarizations and molecular orientations is 
not performed here. Although this would be important 
to describe realistic experimental situations, the averaged 
results generally show smaller discrepancies (not shown). 

At T = 0 fs, the two results show the same profile 
for diagonal and cross peaks and, indeed, the agreement 
is excellent. This shows that the CMRT correctly mod¬ 
els the coherence times and the system-bath correlations 
created during the simulation. With increasing waiting 
time, the same coherent dynamics is found for both the 
diagonal and the cross peaks and even can be inspected 
by eye. However, some disagreement is observed at long 
waiting time T = 500 fs. The diagonal peak B in left 
figure (CMRT-I-PMA) shows a somewhat reduced ampli¬ 
tude as compared to the right Hgure (QUAPI). 

For a more rehned comparison, the amplitudes of the 
diagonal and cross peaks (A, B, C and D) are plotted 
against the waiting time in Fig. [3] and Fig. ID In Fig. [3l 
the population dynamics of the diagonal peaks A (top) 
and B (bottom) calculated by CMRT-I-PMA from 0 to 
1000 fs is shown and compared to the QUAPI result. 
We find that the CMRT-I-PMA provides reasonably ac¬ 
curate results for the population transfer and the oscil¬ 
lation period. However, the amplitude of peak B decays 
slightly faster in the approximate results as compared to 
the QUAPI data. Moreover, both yield different station¬ 
ary states. In addition, the phase of the oscillations is 
slightly shifted. For the comparison of the cross peaks, 
the oscillatory behavior of peaks C and D is plotted in 
Fig- HI Cross peak C shows a similar oscillatory behavior 
but the two approaches yield different stationary states. 
Peak D shows a only slightly shifted phase of the oscil¬ 
latory behavior. Such a phase shift was also observed 
in the population dynamics of the FMO complex shown 
above. The phase shift might be due to the neglect of 
imaginary parts in the Redfield relaxation tensor. 

In order to further assess the reliability of the 
CMRT-I-PMA, we have repeated the calculations for a 
larger reorganization energy, i.e., for A=100 cm“^ (with 
7=100 cm“^ kept unchanged). 2D spectra were again 
calculated by both approaches and the amplitude of the 
labeled peaks were extracted. Their time-dependence 
is plotted in Figs. [S] and IHl CMRT-I-PMA still yields 
quantitative agreement with the QUAPI result except 
for the behavior of the damping. The stronger system- 
bath coupling results in faster damping (diagonal peak 
A) and also in an increased difference between QUAPI 
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FIG. 2. Two-dimensional photon-echo spectra of a dimer 
model calculated by CMRT-I-PMA (left) and QUAPI (right) 
for different waiting times as indicated. The Debye spec¬ 
tral density was used for the calculation with the parameters 
A = 50 cm“^, 7 = 100 cm“^ and the temperature was set to 
T = 77 K. 


and CMRT-I-PMA as compared to the weaker coupling 
with A=50 cm“^. 



T(fs) 

FIG. 3. Population dynamics of the labeled diagonal peaks 
(A, B) extracted from the underlying sequence of 2D maps. 
The two approaches yield the same oscillation period. The 
diagonal peak B obatained from CMRT-I-PMA decays faster 
as compared to the QUAPI result. The oscillation periods can 
be extracted by data htting and are: CMRT-I-PMA; llOfs, 
QUAPI: 99fs. 


V. DISCUSSION - CMRT+PMA VERSUS 
QUAPI 

From the above comparison of the results obtaiired by 
both approaches, we observe that the discrepancies found 
in the 2D calculations are more pronounced than in the 
dynamics of the populations. Put differently, nonlinear 
2D spectra are more sensitive to assess the performance 
and reliability of approximate theoretical approaches. In 
order to understand this, we point out two fundamental 
differences between 2D spectra and the population dy¬ 
namics. First, entanglement between the system and the 
bath leads to initial correlations at the beginning of the 
waiting time window, which are absent in the calcula¬ 
tion of the population dynamics. Second, two-exciton 
states contribute to the 2D spectra during the detec¬ 
tion time, and interference between positive and nega¬ 
tive peaks changes the observed amplitudes. This shows 
that one can not understand the reliability of a method 
to simulate correct 2D spectra by calculating population 
dynamics alone. Our current framework, in which we 
use the combined CMRT-I-PMA and compare the results 
with QUAPI, is well suited to show the performance of 
these methods in understanding 2D spectra directly. 

In more detail, we have observed three noticeable dis¬ 
crepancies of the CMRT-I-PMA as compared to QUAPI: 
i) shifted oscillation phase of peak intensities, ii) a slightly 
faster decay, and, iii) a different amplitude of peaks B and 
C for long waiting times. 

For the explanation of the shifted oscillation observed 



























T(fs) 

FIG. 4. Coherent oscillations of labeled cross peaks (C, 
D) extracted from 2D maps. Cross peak C obained by 
CMRT+PMA shows the same oscillatory behavior, but with 
a somewhat smaller amplitude. 


in 2D simulations of the CMRT+PMA, we need to notice 
that Eqs.[THl[Il]provide the analytic result for a monomer 
(two-level system), and that this has been proven by com¬ 
paring to QUAPI [d^. However, CMRT yields a shifted 
period for the dimer model. The mismatch is mainly 
caused by the population transfer term R(t) since there 
is no population transfer term in the monomer model. In 
this paper, the population transfer rate was calculated 
by the cumulant expansion in Eq. [TT] [l8j | and we only 
took the real part. It is well known that the imaginary 
part dominates the phase of the oscillations [d^. So, 
most likely, the shifted oscillation is mainly caused by 
the real-value approximation of the population transfer 
rates. 

Then, the population transfer term is also derived 
based on the second-order perturbation approximation, 
which is one of the reasons for the explanation of the 
slightly too fast decay of the oscillations found in CMRT 
calculations. Eurthermore, the secular approximation 
was used to separate the population dynamics and the 
dephasing process in Eqs. HSl [121 This also contributes 
to the discrepancy in the decay rate, since it neglects the 
interference between population transfer and coherence 
dephasing. 

A relatively small amplitude of peak B and C was 
found in Eig. [3] and Fig. |d] and it also can be observed 
by eye in the 2D map for the long waiting time T = 500 
fs. We observe that peaks B and C are mainly formed 
by one positive (red) peak and overlap with a negative 
(blue) peak in the 2D spectrum (T = 500 fs). Therefore, 
the amplitude of those peaks mainly depends on the over¬ 
lap of two peaks. In the QUAPI result, the two peaks 
are clearly separated with a larger spectral distance than 



T(fs) 

FIG. 5. Amplitude of the diagonal peaks A and B for 
a stronger system-bath coupling A=100 cm“^ (with 7=100 
cm“^ unchanged). CMRT+PMA calculations yield a faster 
decay (A) as compared to the QUAPI result (decay rate ex¬ 
tracted from a fit: CMRT+PMA: 81 fs, QUAPI: 146 fs). 


in the CMRT result and this leads to the larger ampli¬ 
tude of peaks B and C in the 2D spectrum calculated 
with QUAPI. It indicates that, besides the shifted os¬ 
cillation and the faster decay of the oscillation, CMRT 
does not properly account for the reorganization energy 
by the heat bath (diagonal peaks show slightly different 
positions in the 2D map: —190 cm“^ and 190 cm“^ for 
CMRT and —180 cm“^ and 200 cm“^ for QUAPI). In 
the CMRT, the reorganization energy is included in the 
diagonal part of the Hamiltonian by Eq. |3l where it just 
brings in a shift of the excitonic transition frequency 
by the renormalization term Eq.|3]and does not affect the 
dynamics of the off-diagonal terms in the Hamiltonian in 
Eq. (HI). 

On the basis of a clear physical meaning (population 
transfer and dephasing terms) and for the purpose of 
an efficient and fast calculation, the secular approxima¬ 
tion and the second-order perturbation theory were ap¬ 
plied to construct the CMRT. On the one hand, the sec¬ 
ular approximation leads to a separation of the popu¬ 
lation dynamics and dephasing process and avoids any 
complicated interaction terms between diagonal and off- 
diagonal parts in the equation of motion. On the other 
hand, the second-order perturbation theory simplifies the 
population transfers. It is possible to improve the equa¬ 
tion by including higher orders. However, this renders 
the equation considerably more complicated and requires 
more computational resources for the simulation and it 
is a priori unclear how much this improves the accuracy. 
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T(fs) 

FIG. 6. Oscillations of the cross peaks (C, D) with for a 
stronger system-bath interaction (same parameters as in Fig. 
[5|). The cross-peak C calculated by CMRT-I-PMA yields a 
faster decay and a smaller amplitude (decay rates as obtained 
from a fit: CMRT-I-PMA: 62 fs, QUAPI: 177 fs). 


VI. CONCLUSIONS 

In this paper we present the CMRT and compare it 
in more detail to the QUAPI method by calculating 
the population dynamics of selected FMO exciton sites 
and the 2D-spectrum of a model dimer. We found that 
CMRT provides numerical reliable results as compared 
to numerically exact QUAPI calculations for both the 
population dynamics and 2D spectra, as long as the re¬ 


organization energy is not too large compared to the typ¬ 
ical energy gap of the system. Most importantly, it re¬ 
quires smaller computational efforts and orders of mag¬ 
nitude shorter calculation times. It provides us with an 
efficient approach to study the energy transfer in super- 
large molecular complexes and to perform complicated 
2D simulations. 

We found that the 2D profile calculated from 
CMRT-I-PMA agrees well with the corresponding QUAPI 
results. For a quantitative comparison, the amplitudes of 
diagonal and cross peaks were extracted from 2D maps 
and compared to those calculated by QUAPI. Quantita¬ 
tive agreement was found. We observe some discrepan¬ 
cies. In particular, oscillations are shifted, they decay 
slightly faster, and positions of peaks are slightly shifted. 
This becomes more serious if the reorganization energy 
is increased, which is also the case in Ref. (d^. l45l|. 

The simulation protocol developed here can be used for 
arbitrary forms of the spectral density. One can envision 
an approach where the CMRT method is used to simulate 
super-large complexes, while numerically exact methods 
such as QUAPI play a role in benchmarking the accuracy 
of the simulations of smaller systems. 
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